AIM : To Visualize the satellite products (SST, KPAR) as well as the in situ sites in Lyttelton harbour.
Document Note : The maps, plots and app within this document are interactive so make sure you give them a play like zooming in and out in the maps but also on the plots. Clicking on the legend allows to only select and display the time series needed.
The KPAR dataset consists of monthly means of the attenuation coefficient estimated using the absorption and backscattering coefficient using MODIS-AQUA measurements of the emergent flux, radiance leaving the water. The dataset runs from July 2002 to March 2019.
## Read site sheet
site <- read.csv(file=paste0(getwd(),'/Coords_sites.csv'))
##
## KPAR stack and mean
kpar_mean <- raster('A200207_201903_KPAR_Mean_Lyttelton_QAA.tif')
kpar <- stack('A200207_201903_KPAR_MO_Lyttelton_QAA.tif')
##
## Bathy of the zone NIWA grid
## Download Bathy, plot it with sites and transform as contour
bathy <- raster('Bathy_Lyttelton.tif')
bathy[bathy>0] <- NA
bathycontour <- rasterToContour(bathy)
bathycontour <- spTransform(bathycontour,crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"))
##
## Plot KPAR mean and sites on leaflet map
pal <- colorNumeric(rev(rainbow(10)), values(kpar_mean),
na.color = "transparent")
m <- leaflet(site) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
addTiles() %>%# Print the map
addScaleBar(position = "bottomright",options = scaleBarOptions(imperial=F)) %>%
addMouseCoordinates() %>%
addMarkers(site, lat = ~lat,lng = ~lon, popup = ~name) %>%
addRasterImage(kpar_mean,layerId = "Kd", col=pal, opacity = 0.8,project=T) %>%
addLegend(pal = pal, values = values(kpar_mean),title = "Kd (m-1)") %>%
addImageQuery(x=kpar_mean,layerId = "Kd",type='click',project=T,position='bottomleft') %>%
addPolylines(data=bathycontour,color = "black",opacity=0.2,popup = bathycontour$level)
m
Figure 1 -Mean KPAR product around Lyttelton harbour, NZ.
Key observations:
kpar_NA_sum <- raster('A200207_201903_KPAR_PixAvailability_Lyttelton.tif')
kpar_NA_sum_proj <- projectRaster(kpar_NA_sum,crs=crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"),method='ngb')
pal <- colorNumeric(rev(heat.colors(100)), values(kpar_NA_sum_proj),
na.color = "transparent")
m <- leaflet(site) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
addTiles() %>%# Print the map
addScaleBar(position = "bottomright",options = scaleBarOptions(imperial=F)) %>%
addMouseCoordinates() %>%
addMarkers(site, lat = ~lat,lng = ~lon, popup = ~name) %>%
addRasterImage(kpar_NA_sum_proj,layerId = "Number of pixels available", col=pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(kpar_NA_sum_proj),title = "Number of pixels available") %>%
addImageQuery(x=kpar_NA_sum_proj,layerId = "Number of pixels available",type='click',position='bottomleft') %>%
addPolylines(data=bathycontour,color = "black",opacity=0.2,popup = bathycontour$level)
m
Figure 2 - Availability of pixels for the KPAR product around Lyttelton harbour, NZ.
Key observations:
The SST dataset runs from July 2002 to March 2019 and consists of monthly means. I don’t know from which satellite it is from but could retrieve that at the occasion.
## SST stack and mean
sst <- stack('A200207_201903_SST_MO_Lyttelton.tif')
sst_mean <- raster('A200207_201903_SST_Mean_Lyttelton.tif')
##
## Plot SST mean and sites on leaflet map
pal <- colorNumeric(rev(rainbow(100)), values(sst_mean),
na.color = "transparent")
m <- leaflet(site) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
addTiles() %>%# Print the map
addScaleBar(position = "bottomright",options = scaleBarOptions(imperial=F)) %>%
addMouseCoordinates() %>%
addMarkers(site, lat = ~lat,lng = ~lon, popup = ~name) %>%
addRasterImage(sst_mean,layerId = "SST (degC)", col=pal, opacity = 0.8,project=T) %>%
addLegend(pal = pal, values = values(sst_mean),title = "SST (degC)") %>%
addImageQuery(x=sst_mean,layerId = "SST (degC)",type='click',project=T,position='bottomleft') %>%
addPolylines(data=bathycontour,color = "black",opacity=0.2,popup = bathycontour$level)
m
Figure 3 - Mean SST product around Lyttelton harbour, NZ.
Key observations:
sst_NA_sum <- raster('A200207_201903_SST_PixAvailability_Lyttelton.tif')
pal <- colorNumeric(rev(heat.colors(100)), values(sst_NA_sum),
na.color = "transparent")
m <- leaflet(site) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
addTiles() %>%# Print the map
addScaleBar(position = "bottomright",options = scaleBarOptions(imperial=F)) %>%
addMouseCoordinates() %>%
addMarkers(site, lat = ~lat,lng = ~lon, popup = ~name) %>%
addRasterImage(sst_NA_sum,layerId = "Number of pixels available", col=pal, opacity = 0.8,project=T) %>%
addLegend(pal = pal, values = values(sst_NA_sum),title = "Number of pixels available") %>%
addImageQuery(x=sst_NA_sum,layerId = "Number of pixels available",type='click',position='bottomleft',project=T) %>%
addPolylines(data=bathycontour,color = "black",opacity=0.2,popup = bathycontour$level)
m
Figure 4 - Availability of pixels for the SST product around Lyttelton harbour, NZ.
Key observations:
Different scenarios:
In order to have estimates of KPAR and SST for each sites, we investigate here the scenario of extracting the time series of the pixel where the site is. This approach is supposed to fail as the sites are very coastal so prone to pixels unavailability due to cloud coverage as well as the land adjacency effect which leads to failure in the atmospheric correction process.
## Work on Lat/Lon of site coordinates, reptroject to utm
# Reprojection of sites lon/lat into tmerc (EBED crs)
LatLong_site <- data.frame(Y=site$lat,X=site$lon)
names(LatLong_site) <- c("Y","X")
coordinates(LatLong_site) <- ~ X + Y # longitude first
proj4string(LatLong_site) <- crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0")
utm_sites <- spTransform(LatLong_site, crs( "+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
##
## Extract values of pixel availability at sites location
kpar_NA_sum_velox <- velox('A200207_201903_KPAR_PixAvailability_Lyttelton.tif')
kpar_NA_sum_sites <- kpar_NA_sum_velox$extract_points(utm_sites)
sst_NA_sum_velox <- velox('A200207_201903_SST_PixAvailability_Lyttelton.tif')
sst_NA_sum_sites <- sst_NA_sum_velox$extract_points(utm_sites)
##
## KPAR
kpar_velox <- velox('A200207_201903_KPAR_MO_Lyttelton_QAA.tif')
kpar_sites <- kpar_velox$extract_points(utm_sites)
kpar_sites_df <- data.frame(t(kpar_sites))
colnames(kpar_sites_df) <- site$name
rownames(kpar_sites_df) <- NULL
dateseq <- seq.Date(as.Date("2002/7/1"), by = "month", length.out = dim(kpar_sites)[2])
kpar_sites_df$Date <- dateseq
kpar_sites_tb <- as_tibble(kpar_sites_df)
kpar_sites_tb <- kpar_sites_tb %>% pivot_longer(-Date,names_to="Site",values_to='Kpar') %>% group_by(Date)
q <- ggplot(kpar_sites_tb,aes(Date,Kpar,color=Site)) + geom_point() + geom_line() + labs(y="Kpar (/m)", x = "Date")
ggplotly(q)
Figure 5 - Time series of KPAR (/m) at the 21 sites around Lyttelton harbour, NZ. The values are extracted from the pixels where the sites are (Scenario 1).
Relevant observations:
## KPAR mean, sd, median, min and max values over period
summary_kpar <- kpar_sites_tb %>% group_by(Site) %>% summarise(Mean=mean(Kpar,na.rm=T),Sd=sd(Kpar,na.rm=T),Median=median(Kpar,na.rm=T),Min=min(Kpar,na.rm=T),Max=max(Kpar,na.rm=T))
## `summarise()` ungrouping output (override with `.groups` argument)
summary_kpar$PixAv <- kpar_NA_sum_sites
kable(summary_kpar)
| Site | Mean | Sd | Median | Min | Max | PixAv |
|---|---|---|---|---|---|---|
| BP01 | 0.4807346 | 0.1439982 | 0.4621158 | 0.2263836 | 0.9213209 | 86 |
| BP02 | 0.5245741 | 0.1420317 | 0.5049014 | 0.2587270 | 0.9719090 | 83 |
| BP03 | 0.5477195 | 0.1388396 | 0.5400223 | 0.3038899 | 0.9552073 | 48 |
| BP05 | 0.3739728 | 0.1198098 | 0.3444114 | 0.2331560 | 0.7454379 | 47 |
| BP06 | 0.4466732 | 0.1221743 | 0.4181812 | 0.2577299 | 0.9109438 | 106 |
| BP08 | 0.4077619 | 0.1305797 | 0.3795969 | 0.1555093 | 1.0791025 | 90 |
| BP10 | 0.4126234 | 0.1417000 | 0.3839348 | 0.2120700 | 0.8721717 | 67 |
| BP13 | 0.4603269 | 0.1148118 | 0.4427282 | 0.2227825 | 0.7604557 | 90 |
| BP14 | 0.6887431 | 0.2607180 | 0.6435227 | 0.2820656 | 1.4770980 | 70 |
| LH01 | 0.5052427 | 0.1413214 | 0.5029341 | 0.2142351 | 0.8423380 | 73 |
| LH02 | 0.4657920 | 0.1432052 | 0.4895723 | 0.2473924 | 0.6401853 | 9 |
| LH07 | 0.6919909 | 0.2042824 | 0.6676622 | 0.3572536 | 1.4077549 | 60 |
| LH10 | 0.6233862 | 0.1094229 | 0.6374094 | 0.4760661 | 0.7800066 | 7 |
| PB02 | NaN | NA | NA | Inf | -Inf | 0 |
| PB03 | NaN | NA | NA | Inf | -Inf | 0 |
| PB10 | NaN | NA | NA | Inf | -Inf | 0 |
| PB11 | NaN | NA | NA | Inf | -Inf | 0 |
| PL02 | 0.8797934 | 0.2362799 | 0.9601306 | 0.6019558 | 1.1715964 | 5 |
| PL03 | NaN | NA | NA | Inf | -Inf | 0 |
| PL14 | NaN | NA | NA | Inf | -Inf | 0 |
| PL16 | 0.8977092 | 0.2799484 | 0.9074726 | 0.5237424 | 1.2508827 | 8 |
sst_velox <- velox('A200207_201903_SST_MO_Lyttelton.tif')
sst_sites <- sst_velox$extract_points(utm_sites)
sst_sites_df <- data.frame(t(sst_sites))
colnames(sst_sites_df) <- site$name
rownames(sst_sites_df) <- NULL
dateseq <- seq.Date(as.Date("2002/7/1"), by = "month", length.out = dim(sst_sites)[2])
sst_sites_df$Date <- dateseq
sst_sites_tb <- as_tibble(sst_sites_df)
sst_sites_tb <- sst_sites_tb %>% pivot_longer(-Date,names_to="Site",values_to='SST') %>% group_by(Date)
p <- ggplot(sst_sites_tb,aes(Date,SST,color=Site)) + geom_point() + geom_line() + labs(y="SST (degC)", x = "Date")
ggplotly(p)
Figure 6 - Time series of SST (degC) at the 21 sites around Lyttelton harbour, NZ. The values are extracted from the pixels where the sites are (Scenario 1).
Relevant observations:
## SST mean, sd, median, min and max values over period
summary_sst <- sst_sites_tb %>% group_by(Site) %>% summarise(Mean=mean(SST,na.rm=T),Sd=sd(SST,na.rm=T),Median=median(SST,na.rm=T),Min=min(SST,na.rm=T),Max=max(SST,na.rm=T))
## `summarise()` ungrouping output (override with `.groups` argument)
summary_sst$PixAv <- sst_NA_sum_sites
kable(summary_sst)
| Site | Mean | Sd | Median | Min | Max | PixAv |
|---|---|---|---|---|---|---|
| BP01 | 16.72645 | 3.505492 | 16.85250 | 8.695 | 24.035 | 100 |
| BP02 | 16.93592 | 3.504716 | 17.20750 | 9.870 | 23.340 | 102 |
| BP03 | 16.67018 | 3.449685 | 17.54500 | 8.170 | 23.110 | 51 |
| BP05 | 14.50368 | 3.624410 | 14.19000 | 8.715 | 19.910 | 19 |
| BP06 | 16.46508 | 3.629748 | 17.31250 | 8.520 | 23.965 | 118 |
| BP08 | 15.74805 | 3.667408 | 16.24000 | 7.805 | 22.180 | 107 |
| BP10 | 16.96604 | 3.810822 | 17.72750 | 8.300 | 23.660 | 78 |
| BP13 | 14.59495 | 4.071016 | 14.70000 | 7.100 | 22.625 | 77 |
| BP14 | 16.18449 | 4.400136 | 16.31474 | 9.260 | 25.740 | 50 |
| LH01 | 14.94987 | 4.823432 | 15.33750 | 7.460 | 25.455 | 74 |
| LH02 | 13.81568 | 5.620958 | 12.02000 | 5.640 | 24.265 | 26 |
| LH07 | 18.14572 | 4.455413 | 18.60000 | 8.575 | 26.540 | 79 |
| LH10 | 22.14909 | 3.711519 | 22.14909 | 12.875 | 30.245 | 23 |
| PB02 | NaN | NA | NA | Inf | -Inf | 0 |
| PB03 | NaN | NA | NA | Inf | -Inf | 0 |
| PB10 | NaN | NA | NA | Inf | -Inf | 0 |
| PB11 | NaN | NA | NA | Inf | -Inf | 0 |
| PL02 | 17.39167 | 4.701594 | 16.73083 | 9.650 | 24.750 | 10 |
| PL03 | NaN | NA | NA | Inf | -Inf | 0 |
| PL14 | 7.54500 | 0.000000 | 7.54500 | 7.545 | 7.545 | 2 |
| PL16 | 18.87106 | 4.213544 | 20.63000 | 9.630 | 24.130 | 27 |
In order to have estimates of KPAR and SST for each sites, we investigate here the scenario of extracting the time series from a pixel with the maximum availability (using the pixel availability raster) within a cluster of 16 pixels adjacent to the site location (17 pixels including site, see what’s the structure of the cluster).
## KPAR
ngb_pix_maxPixAv_index_save <- c()
sites_kpar_mean_save <- c()
for (i in 1:dim(site)[1]) {
sites_tmerc <- data.frame(utm_sites)[i,]
sites_kpar_mean <- raster::extract(kpar_mean,sites_tmerc,cellnumbers=T) #Get cell number of pixel where site
ngb_pixels <- raster::adjacent(kpar_mean,cells=sites_kpar_mean[,1],pairs=F,directions=16,include=T)
ngb_pix_maxPixAv_index <- ngb_pixels[which(kpar_NA_sum[ngb_pixels]==max(kpar_NA_sum[ngb_pixels]))] #Return index of pixel in cluster of 16 pixels around that have max pixel availability
# Issues with PB10 and PB3 -> Seems like no pixels available in vicinity
if (length(ngb_pix_maxPixAv_index)>1) {
#Manually solve the issue affecting original index to pixel
ngb_pix_maxPixAv_index <- sites_kpar_mean[,1]
}
ngb_pix_maxPixAv_index_save <- c(ngb_pix_maxPixAv_index_save,ngb_pix_maxPixAv_index)
}
# Position of the "new" pixels where extractions of scenario 2 is from
pix_index_maxPixAv <- xyFromCell(kpar_mean,cell=ngb_pix_maxPixAv_index_save,spatial=F)#Convert to sp feature, to use velox::
Tmerc_site_s2 <- data.frame(Y=pix_index_maxPixAv[,2],X=pix_index_maxPixAv[,1])
coordinates(Tmerc_site_s2) <- ~ X + Y # longitude first
proj4string(Tmerc_site_s2) <- crs("+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
LatLong_sites_s2 <- spTransform(Tmerc_site_s2,crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"))
LatLong_sites_s2$name <- paste0(site$name,"_s2")
LatLong_sites_s2_df <- data.frame(LatLong_sites_s2)
# Extract time series from "new" pixels
pix_index_maxPixAv_sp <- xyFromCell(kpar_mean,cell=ngb_pix_maxPixAv_index_save,spatial=T)#Convert to sp feature, to use velox::
kpar_sites_maxPixAv <- kpar_velox$extract_points(pix_index_maxPixAv_sp)#Classic extraction process
kpar_sites_maxPixAv_df <- data.frame(t(kpar_sites_maxPixAv))
colnames(kpar_sites_maxPixAv_df) <- LatLong_sites_s2$name#paste0(site$name[1],'_maxPixAv')
rownames(kpar_sites_maxPixAv_df) <- NULL
dateseq <- seq.Date(as.Date("2002/7/1"), by = "month", length.out = dim(kpar_sites_maxPixAv_df)[1])
kpar_sites_maxPixAv_df$Date <- dateseq
kpar_sites_maxPixAv_tb <- as_tibble(kpar_sites_maxPixAv_df)
kpar_sites_maxPixAv_tb <- kpar_sites_maxPixAv_tb %>% pivot_longer(-Date,names_to="Site",values_to='Kpar') %>% group_by(Date)
q <- ggplot(kpar_sites_maxPixAv_tb,aes(Date,Kpar,color=Site)) + geom_point() + geom_line() + labs(y="Kpar (/m)", x = "Date")
ggplotly(q)
Figure 7 - Time series of KPAR (/m) at the 21 sites around Lyttelton harbour, NZ. The values are extracted from the pixels with the highest availability in a cluster of 16 pixels adjacent to each sites. Scenario 2.
m <- leaflet(data=LatLong_sites_s2_df) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
addTiles() %>%# Print the map
addMarkers(LatLong_sites_s2_df, lat = ~Y,lng = ~X, popup = ~name)
m
Figure 8 - Position of the Pixels where time series are extracted in the scenario 2 case.
Relevant observations:
## KPAR mean, sd, median, min and max values over period
summary_kpar_s2 <- kpar_sites_maxPixAv_tb %>% group_by(Site) %>% summarise(Mean=mean(Kpar,na.rm=T),Sd=sd(Kpar,na.rm=T),Median=median(Kpar,na.rm=T),Min=min(Kpar,na.rm=T),Max=max(Kpar,na.rm=T))
## `summarise()` ungrouping output (override with `.groups` argument)
summary_kpar_s2$PixAv <- kpar_NA_sum[ngb_pix_maxPixAv_index_save]
kable(summary_kpar_s2)
| Site | Mean | Sd | Median | Min | Max | PixAv |
|---|---|---|---|---|---|---|
| BP01_s2 | 0.3736091 | 0.0757962 | 0.3536769 | 0.2344388 | 0.7235438 | 168 |
| BP02_s2 | 0.3566422 | 0.0690696 | 0.3449551 | 0.2511407 | 0.7337249 | 184 |
| BP03_s2 | 0.3382366 | 0.0698441 | 0.3276764 | 0.2026570 | 0.6407034 | 169 |
| BP05_s2 | 0.3479043 | 0.0756010 | 0.3367585 | 0.2061284 | 0.6099731 | 163 |
| BP06_s2 | 0.3381993 | 0.0727750 | 0.3218081 | 0.2134089 | 0.5809690 | 190 |
| BP08_s2 | 0.3111466 | 0.0691485 | 0.2943282 | 0.1995939 | 0.6331935 | 190 |
| BP10_s2 | 0.3616482 | 0.0944326 | 0.3516572 | 0.1610703 | 0.9613709 | 185 |
| BP13_s2 | 0.3749535 | 0.0785294 | 0.3640774 | 0.1974660 | 0.6571044 | 180 |
| BP14_s2 | 0.3566313 | 0.0928171 | 0.3355087 | 0.1544371 | 0.6814184 | 177 |
| LH01_s2 | 0.4792453 | 0.1089673 | 0.4638782 | 0.2742959 | 0.8294032 | 138 |
| LH02_s2 | 0.5040968 | 0.1202039 | 0.4734725 | 0.2543523 | 0.9888393 | 145 |
| LH07_s2 | 0.5214943 | 0.1296387 | 0.4985805 | 0.2543523 | 0.9773747 | 146 |
| LH10_s2 | 0.4840434 | 0.1156862 | 0.4757362 | 0.2742959 | 0.8221809 | 136 |
| PB02_s2 | 0.4004347 | 0.0580121 | 0.4073745 | 0.3243970 | 0.4625929 | 4 |
| PB03_s2 | NaN | NA | NA | Inf | -Inf | 0 |
| PB10_s2 | NaN | NA | NA | Inf | -Inf | 0 |
| PB11_s2 | 1.1407515 | NA | 1.1407515 | 1.1407515 | 1.1407515 | 1 |
| PL02_s2 | 0.6564987 | 0.2188256 | 0.6817867 | 0.2699000 | 1.0174295 | 41 |
| PL03_s2 | 1.0310847 | 0.2605226 | 1.1485419 | 0.6019558 | 1.3445294 | 11 |
| PL14_s2 | 1.0310847 | 0.2605226 | 1.1485419 | 0.6019558 | 1.3445294 | 11 |
| PL16_s2 | 0.4569534 | 0.1356378 | 0.4164062 | 0.2840720 | 0.9177421 | 87 |
## SST
ngb_pix_maxPixAv_index_save <- c()
sites_sst_mean_save <- c()
for (i in 1:dim(site)[1]) {
sites_tmerc <- data.frame(utm_sites)[i,]
sites_sst_mean <- raster::extract(kpar_mean,sites_tmerc,cellnumbers=T) #Get cell number of pixel where site
ngb_pixels <- raster::adjacent(sst_mean,cells=sites_sst_mean[,1],pairs=F,directions=16,include=T)
ngb_pix_maxPixAv_index <- ngb_pixels[which(sst_NA_sum[ngb_pixels]==max(sst_NA_sum[ngb_pixels]))] #Return index of pixel in cluster of 16 pixels around that have max pixel availability
# Issues with PB10 and PB3 -> Seems like no pixels available in vicinity
if (length(ngb_pix_maxPixAv_index)>1) {
#Manually solve the issue affecting original index to pixel
ngb_pix_maxPixAv_index <- sites_sst_mean[,1]
}
ngb_pix_maxPixAv_index_save <- c(ngb_pix_maxPixAv_index_save,ngb_pix_maxPixAv_index)
}
# Position of the "new" pixels where extractions of scenario 2 is from
pix_index_maxPixAv_sst <- xyFromCell(sst_mean,cell=ngb_pix_maxPixAv_index_save,spatial=F)#Convert to sp feature, to use velox::
Tmerc_site_s2_sst <- data.frame(Y=pix_index_maxPixAv_sst[,2],X=pix_index_maxPixAv_sst[,1])
coordinates(Tmerc_site_s2_sst) <- ~ X + Y # longitude first
proj4string(Tmerc_site_s2_sst) <- crs("+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
LatLong_sites_s2_sst <- spTransform(Tmerc_site_s2_sst,crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"))
LatLong_sites_s2_sst$name <- paste0(site$name,"_s2")
LatLong_sites_s2_sst_df <- data.frame(LatLong_sites_s2_sst)
# Extract time series from "new" pixels
pix_index_maxPixAv_sp <- xyFromCell(sst_mean,cell=ngb_pix_maxPixAv_index_save,spatial=T)#Convert to sp feature, to use velox::
sst_sites_maxPixAv <- sst_velox$extract_points(pix_index_maxPixAv_sp)#Classic extraction process
sst_sites_maxPixAv_df <- data.frame(t(sst_sites_maxPixAv))
colnames(sst_sites_maxPixAv_df) <- LatLong_sites_s2_sst$name#paste0(site$name[1],'_maxPixAv')
rownames(sst_sites_maxPixAv_df) <- NULL
dateseq <- seq.Date(as.Date("2002/7/1"), by = "month", length.out = dim(sst_sites_maxPixAv_df)[1])
sst_sites_maxPixAv_df$Date <- dateseq
sst_sites_maxPixAv_tb <- as_tibble(sst_sites_maxPixAv_df)
sst_sites_maxPixAv_tb <- sst_sites_maxPixAv_tb %>% pivot_longer(-Date,names_to="Site",values_to='SST') %>% group_by(Date)
q <- ggplot(sst_sites_maxPixAv_tb,aes(Date,SST,color=Site)) + geom_point() + geom_line() + labs(y="SST (degC)", x = "Date")
ggplotly(q)
Figure 9 - Time series of SST (degC) at the 21 sites around Lyttelton harbour, NZ. The values are extracted from the pixels with the highest availability in a cluster of 16 pixels adjacent to each sites. Scenario 2.
m <- leaflet(data=LatLong_sites_s2_sst_df) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
addTiles() %>%
addMarkers(LatLong_sites_s2_sst_df, lat = ~Y,lng = ~X, popup = ~name)
m
Figure 10 - Position of the pixels where time series of SST are extracted from in the scenario 2 case.
Relevant observations:
## KPAR mean, sd, median, min and max values over period
# SST mean, sd, median, min and max values over period
summary_sst_s2 <- sst_sites_maxPixAv_tb %>% group_by(Site) %>% summarise(Mean=mean(SST,na.rm=T),Sd=sd(SST,na.rm=T),Median=median(SST,na.rm=T),Min=min(SST,na.rm=T),Max=max(SST,na.rm=T))
## `summarise()` ungrouping output (override with `.groups` argument)
summary_sst_s2$PixAv <- sst_NA_sum[ngb_pix_maxPixAv_index_save]
kable(summary_sst_s2)
| Site | Mean | Sd | Median | Min | Max | PixAv |
|---|---|---|---|---|---|---|
| BP01_s2 | 14.55268 | 3.2726353 | 15.1200 | 7.785000 | 20.520 | 172 |
| BP02_s2 | 16.93592 | 3.5047162 | 17.2075 | 9.870000 | 23.340 | 102 |
| BP03_s2 | 14.35078 | 3.2165821 | 14.6975 | 8.120000 | 20.775 | 166 |
| BP05_s2 | 14.69212 | 3.3178140 | 14.9900 | 8.479999 | 21.245 | 142 |
| BP06_s2 | 13.86971 | 3.3077064 | 13.9600 | 8.264999 | 19.750 | 194 |
| BP08_s2 | 13.94149 | 3.2028660 | 14.2700 | 8.160000 | 19.295 | 187 |
| BP10_s2 | 13.93659 | 3.1511820 | 14.1350 | 8.255000 | 19.590 | 184 |
| BP13_s2 | 14.47864 | 3.5653955 | 14.6300 | 8.264999 | 21.900 | 160 |
| BP14_s2 | 14.44132 | 3.4180549 | 14.6200 | 8.040000 | 21.285 | 157 |
| LH01_s2 | 14.78414 | 3.8983909 | 15.4800 | 7.535000 | 22.210 | 153 |
| LH02_s2 | 14.93382 | 4.3012935 | 15.1625 | 7.610000 | 25.065 | 164 |
| LH07_s2 | 14.82535 | 3.7880783 | 15.3200 | 7.580000 | 22.195 | 155 |
| LH10_s2 | 14.95026 | 3.8885392 | 15.6750 | 7.770000 | 22.660 | 151 |
| PB02_s2 | 13.69000 | 0.3400002 | 13.6900 | 13.349999 | 14.030 | 3 |
| PB03_s2 | NaN | NA | NA | Inf | -Inf | 0 |
| PB10_s2 | NaN | NA | NA | Inf | -Inf | 0 |
| PB11_s2 | 21.83500 | 6.1904860 | 22.7275 | 13.520000 | 28.365 | 4 |
| PL02_s2 | 19.21261 | 4.1895382 | 19.9050 | 9.525000 | 28.535 | 58 |
| PL03_s2 | 19.00000 | 3.9163583 | 18.7625 | 13.145000 | 25.945 | 14 |
| PL14_s2 | 17.37219 | 5.2478802 | 17.8000 | 5.345000 | 25.330 | 17 |
| PL16_s2 | 16.41398 | 3.6061747 | 16.9075 | 8.139999 | 22.940 | 106 |
The following tables show the difference of the values from the scenario 1 minus scenario 2 for each sites. Negative values mean that the scenario 2 leads to an increase in the parameter/value from scenario 1.
## KPAR - Relative difference: s1 - s2
is.na(summary_kpar) <- sapply(summary_kpar, is.infinite)#Replace INF value by NA
is.na(summary_kpar_s2) <- sapply(summary_kpar_s2, is.infinite)#Replace INF value by NA
rdiff_s1s2_kpar <- summary_kpar[,2:7] - summary_kpar_s2[,2:7]
rdiff_s1s2_kpar$Site <- summary_kpar$Site
kable(rdiff_s1s2_kpar)
| Mean | Sd | Median | Min | Max | PixAv | Site |
|---|---|---|---|---|---|---|
| 0.1071255 | 0.0682020 | 0.1084389 | -0.0080553 | 0.1977770 | -82 | BP01 |
| 0.1679319 | 0.0729621 | 0.1599463 | 0.0075863 | 0.2381842 | -101 | BP02 |
| 0.2094828 | 0.0689955 | 0.2123459 | 0.1012329 | 0.3145038 | -121 | BP03 |
| 0.0260685 | 0.0442089 | 0.0076529 | 0.0270276 | 0.1354648 | -116 | BP05 |
| 0.1084739 | 0.0493994 | 0.0963730 | 0.0443210 | 0.3299748 | -84 | BP06 |
| 0.0966152 | 0.0614312 | 0.0852687 | -0.0440846 | 0.4459090 | -100 | BP08 |
| 0.0509752 | 0.0472675 | 0.0322776 | 0.0509997 | -0.0891992 | -118 | BP10 |
| 0.0853734 | 0.0362824 | 0.0786507 | 0.0253166 | 0.1033512 | -90 | BP13 |
| 0.3321118 | 0.1679009 | 0.3080139 | 0.1276285 | 0.7956796 | -107 | BP14 |
| 0.0259974 | 0.0323541 | 0.0390559 | -0.0600609 | 0.0129349 | -65 | LH01 |
| -0.0383048 | 0.0230013 | 0.0160998 | -0.0069599 | -0.3486540 | -136 | LH02 |
| 0.1704966 | 0.0746438 | 0.1690816 | 0.1029013 | 0.4303802 | -86 | LH07 |
| 0.1393428 | -0.0062633 | 0.1616731 | 0.2017702 | -0.0421743 | -129 | LH10 |
| NaN | NA | NA | NA | NA | -4 | PB02 |
| NaN | NA | NA | NA | NA | 0 | PB03 |
| NaN | NA | NA | NA | NA | 0 | PB10 |
| NaN | NA | NA | NA | NA | -1 | PB11 |
| 0.2232947 | 0.0174544 | 0.2783440 | 0.3320558 | 0.1541669 | -36 | PL02 |
| NaN | NA | NA | NA | NA | -11 | PL03 |
| NaN | NA | NA | NA | NA | -11 | PL14 |
| 0.4407557 | 0.1443106 | 0.4910663 | 0.2396705 | 0.3331407 | -79 | PL16 |
##
Key observations:
## SST - Relative difference: s1 - s2
is.na(summary_sst) <- sapply(summary_sst, is.infinite)#Replace INF value by NA
is.na(summary_sst_s2) <- sapply(summary_sst_s2, is.infinite)#Replace INF value by NA
rdiff_s1s2_sst <- summary_sst[,2:7] - summary_sst_s2[,2:7]
rdiff_s1s2_sst$Site <- summary_sst$Site
kable(rdiff_s1s2_sst)
| Mean | Sd | Median | Min | Max | PixAv | Site |
|---|---|---|---|---|---|---|
| 2.1737668 | 0.2328568 | 1.7325006 | 0.9099998 | 3.5149994 | -72 | BP01 |
| 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0 | BP02 |
| 2.3194007 | 0.2331030 | 2.8475003 | 0.0500002 | 2.3349991 | -115 | BP03 |
| -0.1884315 | 0.3065963 | -0.8000002 | 0.2350006 | -1.3349991 | -123 | BP05 |
| 2.5953728 | 0.3220413 | 3.3524995 | 0.2550001 | 4.2150002 | -76 | BP06 |
| 1.8065653 | 0.4645421 | 1.9700003 | -0.3550000 | 2.8850002 | -80 | BP08 |
| 3.0294443 | 0.6596399 | 3.5924993 | 0.0450001 | 4.0699997 | -106 | BP10 |
| 0.1163033 | 0.5056204 | 0.0700002 | -1.1649995 | 0.7250004 | -83 | BP13 |
| 1.7431729 | 0.9820810 | 1.6947432 | 1.2200003 | 4.4549999 | -107 | BP14 |
| 0.1657225 | 0.9250412 | -0.1424999 | -0.0749998 | 3.2450008 | -79 | LH01 |
| -1.1181394 | 1.3196644 | -3.1424999 | -1.9699998 | -0.7999992 | -138 | LH02 |
| 3.3203697 | 0.6673344 | 3.2800007 | 0.9949999 | 4.3449993 | -76 | LH07 |
| 7.1988307 | -0.1770199 | 6.4740906 | 5.1050000 | 7.5849991 | -128 | LH10 |
| NaN | NA | NA | NA | NA | -3 | PB02 |
| NaN | NA | NA | NA | NA | 0 | PB03 |
| NaN | NA | NA | NA | NA | 0 | PB10 |
| NaN | NA | NA | NA | NA | -4 | PB11 |
| -1.8209427 | 0.5120556 | -3.1741657 | 0.1250000 | -3.7849998 | -48 | PL02 |
| NaN | NA | NA | NA | NA | -14 | PL03 |
| -9.8271873 | -5.2478802 | -10.2549996 | 2.1999998 | -17.7850003 | -15 | PL14 |
| 2.4570789 | 0.6073695 | 3.7224998 | 1.4900007 | 1.1900005 | -79 | PL16 |
Key observations: